Introduction to Ambisonics

Fundamentals

Chris Hold, Virtual Acoustics (Aalto University)

Structure

• Motivation
• Discovering Spherical Harmonics
• Ambisonics Encoder and Decoder

From Recording to Playback

Audio Description Formats

Channel Based	Object Based	Scene Based

Scene Based Workflow

Benefits of Scene-based

• Separation of recording, transmission / storage, and playback
• Flexible with multiple options for each stage
• Does not scale with the number of sources

Ambisonics

• Implementation of a scene based format
• "Long" history and connection to other sciences
• Use spherical harmonics as basis functions

Example

Error in callback <function flush_figures at 0x15c95d1c0> (for post_execute), with arguments args (),kwargs {}:

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib_inline/backend_inline.py:126, in flush_figures()
    123 if InlineBackend.instance().close_figures:
    124     # ignore the tracking, just draw and close all figures
    125     try:
--> 126         return show(True)
    127     except Exception as e:
    128         # safely show traceback if in IPython, else raise
    129         ip = get_ipython()

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib_inline/backend_inline.py:90, in show(close, block)
     88 try:
     89     for figure_manager in Gcf.get_all_fig_managers():
---> 90         display(
     91             figure_manager.canvas.figure,
     92             metadata=_fetch_figure_metadata(figure_manager.canvas.figure)
     93         )
     94 finally:
     95     show._to_draw = []

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/IPython/core/display_functions.py:298, in display(include, exclude, metadata, transient, display_id, raw, clear, *objs, **kwargs)
    296     publish_display_data(data=obj, metadata=metadata, **kwargs)
    297 else:
--> 298     format_dict, md_dict = format(obj, include=include, exclude=exclude)
    299     if not format_dict:
    300         # nothing to display (e.g. _ipython_display_ took over)
    301         continue

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/IPython/core/formatters.py:179, in DisplayFormatter.format(self, obj, include, exclude)
    177 md = None
    178 try:
--> 179     data = formatter(obj)
    180 except:
    181     # FIXME: log the exception
    182     raise

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/decorator.py:232, in decorate.<locals>.fun(*args, **kw)
    230 if not kwsyntax:
    231     args, kw = fix(args, kw, sig)
--> 232 return caller(func, *(extras + args), **kw)

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/IPython/core/formatters.py:223, in catch_format_error(method, self, *args, **kwargs)
    221 """show traceback on failed format call"""
    222 try:
--> 223     r = method(self, *args, **kwargs)
    224 except NotImplementedError:
    225     # don't warn on NotImplementedErrors
    226     return self._check_return(None, args[0])

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/IPython/core/formatters.py:340, in BaseFormatter.__call__(self, obj)
    338     pass
    339 else:
--> 340     return printer(obj)
    341 # Finally look for special method names
    342 method = get_real_method(obj, self.print_method)

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/IPython/core/pylabtools.py:152, in print_figure(fig, fmt, bbox_inches, base64, **kwargs)
    149     from matplotlib.backend_bases import FigureCanvasBase
    150     FigureCanvasBase(fig)
--> 152 fig.canvas.print_figure(bytes_io, **kw)
    153 data = bytes_io.getvalue()
    154 if fmt == 'svg':

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/backend_bases.py:2187, in FigureCanvasBase.print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
   2183 try:
   2184     # _get_renderer may change the figure dpi (as vector formats
   2185     # force the figure dpi to 72), so we need to set it again here.
   2186     with cbook._setattr_cm(self.figure, dpi=dpi):
-> 2187         result = print_method(
   2188             filename,
   2189             facecolor=facecolor,
   2190             edgecolor=edgecolor,
   2191             orientation=orientation,
   2192             bbox_inches_restore=_bbox_inches_restore,
   2193             **kwargs)
   2194 finally:
   2195     if bbox_inches and restore_bbox:

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/backend_bases.py:2043, in FigureCanvasBase._switch_canvas_and_return_print_method.<locals>.<lambda>(*args, **kwargs)
   2039     optional_kws = {  # Passed by print_figure for other renderers.
   2040         "dpi", "facecolor", "edgecolor", "orientation",
   2041         "bbox_inches_restore"}
   2042     skip = optional_kws - {*inspect.signature(meth).parameters}
-> 2043     print_method = functools.wraps(meth)(lambda *args, **kwargs: meth(
   2044         *args, **{k: v for k, v in kwargs.items() if k not in skip}))
   2045 else:  # Let third-parties do as they see fit.
   2046     print_method = meth

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/backends/backend_svg.py:1339, in FigureCanvasSVG.print_svg(self, filename, bbox_inches_restore, metadata)
   1334 w, h = width * 72, height * 72
   1335 renderer = MixedModeRenderer(
   1336     self.figure, width, height, dpi,
   1337     RendererSVG(w, h, fh, image_dpi=dpi, metadata=metadata),
   1338     bbox_inches_restore=bbox_inches_restore)
-> 1339 self.figure.draw(renderer)
   1340 renderer.finalize()

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/artist.py:95, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
     93 @wraps(draw)
     94 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 95     result = draw(artist, renderer, *args, **kwargs)
     96     if renderer._rasterizing:
     97         renderer.stop_rasterizing()

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/artist.py:72, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     69     if artist.get_agg_filter() is not None:
     70         renderer.start_filter()
---> 72     return draw(artist, renderer)
     73 finally:
     74     if artist.get_agg_filter() is not None:

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/figure.py:3154, in Figure.draw(self, renderer)
   3151         # ValueError can occur when resizing a window.
   3153 self.patch.draw(renderer)
-> 3154 mimage._draw_list_compositing_images(
   3155     renderer, self, artists, self.suppressComposite)
   3157 for sfig in self.subfigs:
   3158     sfig.draw(renderer)

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/image.py:132, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
    130 if not_composite or not has_images:
    131     for a in artists:
--> 132         a.draw(renderer)
    133 else:
    134     # Composite any adjacent images together
    135     image_group = []

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/artist.py:39, in _prevent_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
     36     renderer.stop_rasterizing()
     37     renderer._rasterizing = False
---> 39 return draw(artist, renderer, *args, **kwargs)

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/projections/polar.py:1036, in PolarAxes.draw(self, renderer)
   1033     self.yaxis.reset_ticks()
   1034     self.yaxis.set_clip_path(self.patch)
-> 1036 super().draw(renderer)

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/artist.py:72, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     69     if artist.get_agg_filter() is not None:
     70         renderer.start_filter()
---> 72     return draw(artist, renderer)
     73 finally:
     74     if artist.get_agg_filter() is not None:

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/axes/_base.py:3034, in _AxesBase.draw(self, renderer)
   3031     for spine in self.spines.values():
   3032         artists.remove(spine)
-> 3034 self._update_title_position(renderer)
   3036 if not self.axison:
   3037     for _axis in self._axis_map.values():

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/axes/_base.py:2978, in _AxesBase._update_title_position(self, renderer)
   2976 top = max(top, bb.ymax)
   2977 if title.get_text():
-> 2978     ax.yaxis.get_tightbbox(renderer)  # update offsetText
   2979     if ax.yaxis.offsetText.get_text():
   2980         bb = ax.yaxis.offsetText.get_tightbbox(renderer)

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/axis.py:1336, in Axis.get_tightbbox(self, renderer, for_layout_only)
   1333     renderer = self.figure._get_renderer()
   1334 ticks_to_draw = self._update_ticks()
-> 1336 self._update_label_position(renderer)
   1338 # go back to just this axis's tick labels
   1339 tlb1, tlb2 = self._get_ticklabel_bboxes(ticks_to_draw, renderer)

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/axis.py:2609, in YAxis._update_label_position(self, renderer)
   2605     return
   2607 # get bounding boxes for this axis and any siblings
   2608 # that have been set by `fig.align_ylabels()`
-> 2609 bboxes, bboxes2 = self._get_tick_boxes_siblings(renderer=renderer)
   2610 x, y = self.label.get_position()
   2611 if self.label_position == 'left':

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/axis.py:2161, in Axis._get_tick_boxes_siblings(self, renderer)
   2159 axis = ax._axis_map[name]
   2160 ticks_to_draw = axis._update_ticks()
-> 2161 tlb, tlb2 = axis._get_ticklabel_bboxes(ticks_to_draw, renderer)
   2162 bboxes.extend(tlb)
   2163 bboxes2.extend(tlb2)

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/axis.py:1315, in Axis._get_ticklabel_bboxes(self, ticks, renderer)
   1313 if renderer is None:
   1314     renderer = self.figure._get_renderer()
-> 1315 return ([tick.label1.get_window_extent(renderer)
   1316          for tick in ticks if tick.label1.get_visible()],
   1317         [tick.label2.get_window_extent(renderer)
   1318          for tick in ticks if tick.label2.get_visible()])

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/axis.py:1315, in <listcomp>(.0)
   1313 if renderer is None:
   1314     renderer = self.figure._get_renderer()
-> 1315 return ([tick.label1.get_window_extent(renderer)
   1316          for tick in ticks if tick.label1.get_visible()],
   1317         [tick.label2.get_window_extent(renderer)
   1318          for tick in ticks if tick.label2.get_visible()])

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/text.py:958, in Text.get_window_extent(self, renderer, dpi)
    956 bbox, info, descent = self._get_layout(self._renderer)
    957 x, y = self.get_unitless_position()
--> 958 x, y = self.get_transform().transform((x, y))
    959 bbox = bbox.translated(x, y)
    960 return bbox

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/transforms.py:1495, in Transform.transform(self, values)
   1492 values = values.reshape((-1, self.input_dims))
   1494 # Transform the values
-> 1495 res = self.transform_affine(self.transform_non_affine(values))
   1497 # Convert the result back to the shape of the input values.
   1498 if ndim == 0:

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/_api/deprecation.py:297, in rename_parameter.<locals>.wrapper(*args, **kwargs)
    292     warn_deprecated(
    293         since, message=f"The {old!r} parameter of {func.__name__}() "
    294         f"has been renamed {new!r} since Matplotlib {since}; support "
    295         f"for the old name will be dropped %(removal)s.")
    296     kwargs[new] = kwargs.pop(old)
--> 297 return func(*args, **kwargs)

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/transforms.py:2419, in CompositeGenericTransform.transform_non_affine(self, values)
   2417     return self._a.transform_non_affine(values)
   2418 else:
-> 2419     return self._b.transform_non_affine(self._a.transform(values))

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/_api/deprecation.py:297, in rename_parameter.<locals>.wrapper(*args, **kwargs)
    292     warn_deprecated(
    293         since, message=f"The {old!r} parameter of {func.__name__}() "
    294         f"has been renamed {new!r} since Matplotlib {since}; support "
    295         f"for the old name will be dropped %(removal)s.")
    296     kwargs[new] = kwargs.pop(old)
--> 297 return func(*args, **kwargs)

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/transforms.py:2417, in CompositeGenericTransform.transform_non_affine(self, values)
   2415     return values
   2416 elif not self._a.is_affine and self._b.is_affine:
-> 2417     return self._a.transform_non_affine(values)
   2418 else:
   2419     return self._b.transform_non_affine(self._a.transform(values))

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/_api/deprecation.py:297, in rename_parameter.<locals>.wrapper(*args, **kwargs)
    292     warn_deprecated(
    293         since, message=f"The {old!r} parameter of {func.__name__}() "
    294         f"has been renamed {new!r} since Matplotlib {since}; support "
    295         f"for the old name will be dropped %(removal)s.")
    296     kwargs[new] = kwargs.pop(old)
--> 297 return func(*args, **kwargs)

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/transforms.py:2419, in CompositeGenericTransform.transform_non_affine(self, values)
   2417     return self._a.transform_non_affine(values)
   2418 else:
-> 2419     return self._b.transform_non_affine(self._a.transform(values))

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/_api/deprecation.py:297, in rename_parameter.<locals>.wrapper(*args, **kwargs)
    292     warn_deprecated(
    293         since, message=f"The {old!r} parameter of {func.__name__}() "
    294         f"has been renamed {new!r} since Matplotlib {since}; support "
    295         f"for the old name will be dropped %(removal)s.")
    296     kwargs[new] = kwargs.pop(old)
--> 297 return func(*args, **kwargs)

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/projections/polar.py:78, in PolarTransform.transform_non_affine(self, values)
     76     r = (r - self._get_rorigin()) * self._axis.get_rsign()
     77 r = np.where(r >= 0, r, np.nan)
---> 78 return np.column_stack([r * np.cos(theta), r * np.sin(theta)])

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/numpy/lib/shape_base.py:652, in column_stack(tup)
    650         arr = array(arr, copy=False, subok=True, ndmin=2).T
    651     arrays.append(arr)
--> 652 return _nx.concatenate(arrays, 1)

KeyboardInterrupt: 

From the Wave Equation to Spherical Harmonics

Wave equation \begin{equation} \nabla^2 p - \frac{1}{c^2} \frac{\partial^2 p}{\partial t^2} = 0 \quad. \end{equation}

In frequency domain, with wave-number $k = \frac{\omega}{c}$ results in the Helmholtz equation \begin{equation} (\nabla^2 + k^2) p = 0 \quad. \end{equation}

Multiple solutions, e.g. a mono chromatic plane wave with amplitude $\hat{A}(\omega)$ in Cartesian coordinates \begin{equation} p(t, x, y, z) = \hat{A}(\omega) e^{i(k_x x + k_y y + k_z z - \omega t)} \quad. \end{equation}

Any solution to the Helmholtz equation can also be expressed in spherical coordinates \begin{equation} p(r, \theta, \phi, k) = \sum_{n = 0}^{\infty} \sum_{m=-n}^{n} \color{olive}{(A_{mn} j_n(kr) + B_{mn} y_n(kr))} \; \color{brown}{Y_{n}^{m}(\theta,\phi)} \quad, \end{equation}

\begin{equation} p(r, \theta, \phi, k) = \sum_{n = 0}^{\infty} \sum_{m=-n}^{n} \color{olive}{C_{mn}(kr)} \; \color{brown}{Y_{n}^{m}(\theta,\phi)} \quad, \end{equation}

With two separable parts:

• radial component : linear combination of spherical Bessel functions of first ($j_n$) and second kind $y_n$
• angular component : spherical harmonics $\color{brown}{Y_{n}^{m}(\theta,\phi)}$

sound field defined on a sphere is fully captured by its spherical harmonics coefficients

E.g. a unit plane wave \begin{equation} p(r, \theta, \phi, k) = \sum_{n = 0}^{\infty} \sum_{m=-n}^{n} 4 \pi i^n j_n(kr) \left[ Y_n^m(\theta_k,\phi_k) \right] ^* Y_{n}^ {m}(\theta,\phi) \quad. \end{equation}

Discovering Spherical Harmonics

\begin{equation} Y_n^m(\theta,\phi)=\color{orange}{\sqrt{\frac{2n+1}{4\pi}\frac{(n-m)!}{(n+m)!}}}\, \color{teal}{P_n^m(\cos(\theta))}\, \color{purple}{e^{im\phi}} \quad, \end{equation}\begin{equation} Y_n^m(\theta,\phi)=\color{orange}{D_{nm}}\, \color{teal}{P_n^m(\cos(\theta))}\, \color{purple}{e^{im\phi}} \quad, \end{equation}

where

• azimuth as $\phi$ and zenith/colatitude as $\theta$.
• $P_n^m$ is the associated Legendre polynomial of order $n$ and degree $m$.

Combined with appropriate scaling give real Spherical Harmonics $Y_{n,m}(\theta,\phi)$ as

\begin{equation} Y_{n,m}(\theta,\phi) = \sqrt{ \frac{(2n+1)}{4\pi} \frac{(n-|m|)!}{(n+|m|)!} } P_{n,|m|}(\cos(\theta)) \begin{cases} \sqrt2\sin(|m|\phi) & \mathrm{if\hspace{0.5em}} m < 0 \quad,\\ 1 & \mathrm{if\hspace{0.5em}} m = 0 \quad,\\ \sqrt2\cos(|m|\phi) & \mathrm{if\hspace{0.5em}} m > 0 \quad. \end{cases} \end{equation}

• Always check: Condon–Shortley phase convention $(-1)^m$

Error in callback <function _draw_all_if_interactive at 0x11220e520> (for post_execute), with arguments args (),kwargs {}:

KeyboardInterrupt

Let's look at the azimuthal component $e^{im\phi}$ and the zenithal component $P_n^m(\cos\theta)$

azimuthal component

zenithal component

Orthonormality

Two functions $\color{violet}f, \color{purple}g$ over a domain $\gamma$ are orthogonal if

\begin{equation} \int_\gamma \color{violet}{f^*(\gamma)} \color{purple}{g(\gamma)} \,\mathrm{d}\gamma = \langle \color{violet}f, \color{purple}g \rangle = 0 \quad,\,\mathrm{for}~f \neq g. \end{equation}

They are also orthonormal if \begin{equation} \color{violet}{\int_\gamma f^*(\gamma) f(\gamma) \,\mathrm{d}\gamma = \int_\gamma|f(\gamma)|^2 \,\mathrm{d}\gamma = \langle f, f \rangle } = 1 \quad, \\ \color{purple}{\int_\gamma g^*(\gamma) g(\gamma) \,\mathrm{d}\gamma = \int_\gamma |g(\gamma)|^2 \,\mathrm{d}\gamma = \langle g, g \rangle} = 1 \quad. \end{equation}

For the Spherical Harmonics:

• the the azimuthal component $e^{im\phi}$ along $\phi$ is orthogonal w.r.t. the degree $m$
• the zenithal component $P_n^m(\cos\theta)$ along $\theta$ is orthogonal w.r.t. the order $n$.

Their product is still orthogonal, and the scaling $\color{orange}{D_{nm}}$ ensures orthonormality such that \begin{equation} \int_\Omega Y_n^m(\Omega)^* \, Y_{n'}^{m'}(\Omega) \,\mathrm{d}\Omega = \langle Y_n^m(\Omega) , Y_{n'}^{m'}(\Omega) \rangle = \delta_{nn'}\delta_{mm'} \quad , \end{equation}

and

\begin{equation} \int_{{\Omega} \in \mathbb{S}^2} |Y_n^m({\Omega})|^2 \mathrm{d}{\Omega} = 1 \quad . \end{equation}

Example

show normality for $Y_0^0$: $$ Y_0^0(\theta,\phi)=\sqrt{\frac{0+1}{4\pi}\frac{(0)!}{(0)!}} P_0^0(\cos(\theta)) e^{i0\phi} = \sqrt{\frac{1}{4\pi}} \quad, $$ hence $$ \int_{{\Omega} \in \mathbb{S}^2} Y_0^0(\theta,\phi)^* Y_0^0(\theta,\phi) \,\mathrm{d}{\Omega} = \int_{{\Omega} \in \mathbb{S}^2} \sqrt{\frac{1}{4\pi}}\sqrt{\frac{1}{4\pi}} \, \mathrm{d}{\Omega} = 4\pi \frac{1}{4\pi} = 1 \quad. $$

Example

Test orthogonality of functions $\color{violet}f, \color{purple}g$

$$ \langle \color{violet}f, \color{purple}g \rangle \overset{?}{=} 0 \quad,\,\mathrm{for}~f \neq g. $$

Test for orthogonality in azimuth <f_azi, g_azi> =  -0.0
Test for orthogonality in zenith <f_zen, g_zen> =  0.0

Spherical Harmonic Transform (SHT)

• Find a scene based encoding for sound field
• We showed that the spherical harmonics $Y_n^m({\Omega})$ are a set of suitable basis functions on the sphere
• We also showed that a sound field (on the sphere) $s({\Omega})$ is fully captured by its spherical harmonics coefficients $\sigma_{nm}$

This can be expressed with $\Omega = [\phi, \theta]$ as the inverse Spherical Harmonic Transform (iSHT) \begin{equation} s({\Omega}) = \sum_{n = 0}^{N=\infty} \sum_{m=-n}^{+n} \sigma_{nm} Y_n^m({\Omega}) \quad. \end{equation}

Spherical harmonics coefficients $\sigma_{nm}$ can be derived with the Spherical Harmonic Transform (SHT) \begin{equation} \sigma_{nm} = \int_{{\Omega} \in \mathbb{S}^2} s({\Omega}) [Y_n^m({\Omega})]^* \mathrm{d}{\Omega} = \langle [Y_n^m({\Omega})] , s({\Omega}) \rangle \quad. \end{equation}

• SHT also refered to as spherical Fourier transform
• Note the transform from a spatially discrete to continuous domain

Spherical Grids

The SHT evaluates the continuous integral over $\Omega$ \begin{equation} \sigma_{nm} = \int_{{\Omega} \in \mathbb{S}^2} s({\Omega}) [Y_n^m({\Omega})]^* \mathrm{d}{\Omega} \quad . \end{equation} Quadrature methods allow evaluation by spherical sampling at certain (weighted) grid points such that \begin{equation} \sigma_{nm} \approx \sum_{q=1}^{Q} w_q s({\Omega}_q) [Y_n^m({\Omega}_q)]^* \quad. \end{equation}

Certain grids with sampling points ${\Omega}_q$ and associated sampling weights $w_q$ have certain properties:

• quadrature grids allow numerical integration of spherical polynomials
• Spherical sampling dictates maximum order
• easy to evaluate are uniform/regular grids, with constant $w_q = \frac{4\pi}{Q}$
• An example are so-called spherical t-designs(t), which allow a SHT up to order $N = \lfloor t/2 \rfloor$

Spherical Dirac

We can show that spherical harmonics are orthogonal (even orthonormal) with \begin{equation} \int_\Omega Y_n^m(\Omega) \, Y_{n'}^{m'}(\Omega) \,\mathrm{d}\Omega = \delta_{nn'}\delta_{mm'} \quad . \end{equation}

Because of their completeness, we can also directly formulate a Dirac function on the sphere as \begin{equation} \sum_{n=0}^{N=\infty} \sum_{m=-n}^n [Y_n^m({\Omega'})]^* Y_n^m(\Omega) = \delta(\Omega - \Omega') \quad, \end{equation}

and therefore the spherical Fourier coefficients $\sigma_{nm}$ \begin{equation} SHT\{\delta(\Omega - \Omega')\} = \int_{{\Omega} \in \mathbb{S}^2} \delta(\Omega - \Omega') \, [Y_n^m({\Omega})]^* \mathrm{d}{\Omega} = [Y_n^m({\Omega'})]^* \quad . \end{equation}

Order-Limitation of Spherical Dirac Pulse

Example

Integrate (order-limited) Dirac $\delta_N(\Omega - \Omega')$ over sphere

$$ \int_{{\Omega} \in \mathbb{S}^2} \delta_N(\Omega - \Omega') \mathrm{d}{\Omega} = \int_{{\Omega} \in \mathbb{S}^2} \color{blue}{\sum_{n=0}^{N} \sum_{m=-n}^n\, [Y_n^m({\Omega'})]^* Y_n^m(\Omega)} \,\mathrm{d}{\Omega} \quad, $$

by discretization with sufficient t-design $$ \int_{{\Omega} \in \mathbb{S}^2} \color{blue}{\sum_{n=0}^{N} \sum_{m=-n}^n\, [Y_n^m({\Omega'})]^* Y_n^m(\Omega)} \,\mathrm{d}{\Omega} = \frac{4\pi}{Q}\sum_{q=1}^{Q} \color{orange}{\sum_{n=0}^{N} \sum_{m=-n}^n\, [Y_n^m({\Omega'})]^* Y_n^m(\Omega_q)} $$

Matrix Notations

Stack the spherical harmonics evaluated at $\Omega$ up to spherical order $N$ as

$$ \mathbf{Y} = \left[ \begin{array}{ccccc} Y_0^0(\Omega[0]) & Y_1^{-1}(\Omega[0]) & Y_1^0(\Omega[0]) & \dots & Y_N^N(\Omega[0]) \\ Y_0^0(\Omega[1]) & Y_1^{-1}(\Omega[1]) & Y_1^0(\Omega[1]) & \dots & Y_N^N(\Omega[1]) \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ Y_0^0(\Omega[Q-1]) & Y_1^{-1}(\Omega[Q-1]) & Y_1^0(\Omega[Q-1]) & \dots & Y_N^N(\Omega[Q-1]) \end{array} \right] $$

such that $\mathbf{Y} : [Q, (N+1)^2]$

• soundfield pressure is real signal -> only real spherical harmonics needed
• orthonormal scaling introduced earlier called N3D by convention (in contrast to SN3D)
• stacking them as above, where $idx_{n,m} = n^2+n+m$, is called ACN

• We can stack a time signal $s(t_0, t_1, \ldots, t_T)$ as vector $\mathbf{s}$.

\begin{equation*} \mathbf{s} = \begin{bmatrix} s(t_0) \\ s(t_1) \\ \vdots \\ s(t_T) \end{bmatrix} \end{equation*}

• This leads to a matrix notation of multiple discrete signals $\mathbf{S} : [T, Q]$.

\begin{equation*} \mathbf{S} = \begin{bmatrix} s_0(t_0) & s_1(t_0) & \cdots & s_Q(t_0)\\ s_0(t_1) & s_1(t_1) & \cdots & s_Q(t_1)\\ \vdots & \vdots& \ddots & \vdots \\ s_0(t_T) & s_1(t_T) & \cdots & s_Q(t_T) \end{bmatrix} \end{equation*}

• Similarly, stacking Ambisonics coefficients $\sigma_n^m(t_0, t_1, \ldots, t_T)$ into $\mathbf{\chi} : [T, (N+1)^2]$

\begin{equation*} \mathbf{\chi} = \begin{bmatrix} \sigma_0^0(t_0) & \sigma_1^{-1}(t_0) & \cdots & \sigma_N^N(t_0) \\ \sigma_0^0(t_1) & \sigma_1^{-1}(t_1) & \cdots & \sigma_N^N(t_1) \\ \vdots & \vdots & \ddots & \vdots \\ \sigma_0^0(t_T) & \sigma_1^{-1}(t_T) & \cdots & \sigma_N^N(t_T) \end{bmatrix} \end{equation*}

We obtain ambisonic signals matrix $\mathbf{\chi} : [T, (N+1)^2]$ from signals $\mathbf{S} : [T, Q]$ by SHT in matrix notation as \begin{equation} \mathbf{\chi}(t) = \mathbf{S}(t) \, \mathrm{diag}(w_q) \, \mathbf{Y} \quad. \end{equation} With the SH basis functions evaluated at $\Omega_q$ as $\mathbf{Y} : [Q, (N+1)^2]$.

Obtaining the discrete signals $s_q(t)$ is a linear combination of SH basis functions evaluated at $\Omega_q$ as $\mathbf{Y} : [Q, (N+1)^2]$ and the SH coefficients. This inverse SHT in matrix notation with ambisonic signals matrix $\mathbf{\chi} : [T, (N+1)^2]$ is

\begin{equation} s_q(t) = \mathbf{\chi}(t) \, \mathbf{Y}^T \quad . \end{equation}

From the Ambisonics Encoder, to Directional Weighting, to a Decoder

Encoder

A single plane wave encoded in direction $\Omega_1$ with signal $\mathbf{s}$ is directly the outer product with the Dirac SH coefficients

\begin{equation} \mathbf{\chi}_\mathrm{PW(\Omega_1)} = \mathbf{s} \, \mathbf{Y}(\Omega_1) \quad. \end{equation}

For multiple sources $Q$, we stack and sum \begin{equation} \mathbf{\chi}_\mathrm{PW(\Omega_Q)} = \sum_{q=1}^Q\mathbf{s}_q \, \mathbf{Y}(\Omega_q) = \mathbf{S} \, \mathbf{Y} \quad. \end{equation}

Directional Weighting

• Directionally filtering a soundfield in direction $\Omega_k$ by weighting
• Weighting in SH domain $w_{nm}$ is an elegant way of beamforming

The simplest beamformer is a spherical Dirac in direction $\Omega_k$ normalized by its energy, i.e. $\mathrm{max}DI$ \begin{equation} w_{nm, \mathrm{maxDI}}(\Omega_k) = \frac{4\pi}{(N+1)^2} Y_{n,m}(\Omega_k) \end{equation}

Other patterns can be achieved by weighting the spherical Fourier spectrum. Axis-symmetric patterns reduce to only a modal weighting $c_n$, such that \begin{equation} w_{nm}(\Omega_k) = c_{n} \, Y_{n,m}(\Omega_k) \end{equation}

E.g. $\mathrm{max}\vec{r}_E$ weights each order $n$ with \begin{equation} c_{n,\,\mathrm{max}\vec{r}_E} = P_n[\cos(\frac{137.9^\circ}{N+1.51})] \quad, \end{equation} with the Legendre polynomials $P_n$ of order $n$.

Or we can define a (SH) Butterworth filter with \begin{equation} c_{n,\,\mathrm{Butterworth}} = \frac{1}{\sqrt{1+(n/n_c)^{2k}}} \quad, \end{equation} with the filter-order $k$ and cut-on SH-order $n_c$.

Decoder

• Beamformer extracts a portion of the soundfield in steering direction $\Omega$
• Pointing a set of beamformers in directions $\Omega_k$ results in a set of signal components
• in case of $\mathrm{max}DI$ proportional to iSHT in $\Omega_k$

\begin{equation} s({\Omega_k}) = \sum_{n = 0}^{N} \sum_{m=-n}^{+n} w_{nm}({\Omega_k}) \, \sigma_{nm} \quad, \end{equation}

or in matrix notation with $\mathbf{S}$ and beamforming weights $\mathbf{c}_n$ and $\mathbf{Y}$ \begin{equation} \mathbf{S} = \mathbf{\chi} \, \mathrm{diag_N}(\mathbf{c}_n) \, \mathbf{Y}^T \quad . \end{equation}

Example

SHD signal of 3 sine signal PWs plus noise -> three beamformers, two pointing to sources

Example

Decode on a t-design(6) (sufficient up to $N = 3$):

• a 3rd order ambisonic signal

Example

Decode on a t-design(6) (sufficient up to $N = 3$):

• a 5th order ambisonic signal
• an 8th order ambisonic signal

Loudspeaker Decoders

# Loudspeaker Setup
ls_dirs = np.array([[-135, -80, -45, 0, 45, 80, 135, -135, -60, -30, 30, 60, 135],
                    [0, 0, 0, 0, 0, 0, 0, 60, 60, 60, 60, 60, 60]])
ls_x, ls_y, ls_z = spa.utils.sph2cart(spa.utils.deg2rad(ls_dirs[0, :]),
                                      spa.utils.deg2rad(90 - ls_dirs[1, :]))
ls_setup = spa.decoder.LoudspeakerSetup(ls_x, ls_y, ls_z)
ls_setup.pop_triangles(normal_limit=85, aperture_limit=90,
                       opening_limit=150)
ls_setup.show()

Face not pointing towards listener: [0 1 2]
Face not pointing towards listener: [0 6 2]
Face not pointing towards listener: [2 5 6]
Face not pointing towards listener: [2 3 4]
Face not pointing towards listener: [2 5 4]

spa.plot.decoder_performance(ls_setup, 'NLS')

spa.plot.decoder_performance(ls_setup, 'VBAP')

spa.plot.decoder_performance(ls_setup, 'SAD')

spa.plot.decoder_performance(ls_setup, 'MAD')

spa.plot.decoder_performance(ls_setup, 'EPAD', N_sph=2)

ls_setup.ambisonics_setup(update_hull=True)
spa.plot.decoder_performance(ls_setup, 'ALLRAD')

ls_setup.ambisonics_setup(update_hull=True)
spa.plot.decoder_performance(ls_setup, 'ALLRAD2')

Binaural Decoding

\begin{equation} s^{l,r}(t) = x(t) * h_{\mathrm{HRIR}}^{l,r}({\Omega}, t) \quad, \end{equation}

where $(*)$ denotes the time-domain convolution operation.

Transforming to the time-frequency domain through the time-domain Fourier transform, further assuming plane-wave components $\bar X (\Omega)$, the ear input signals are given as \begin{equation} S^{l,r}(\omega) = \int_{\Omega} \bar X (\Omega, \omega) H^{l,r}(\Omega, \omega) \,\mathrm{d}\Omega \quad. \end{equation}

\begin{equation} S^{l,r}(\omega) = \sum_{n = 0}^{N} \sum_{m = -n}^{+n} \chi_{nm}(\omega) \breve H_{nm}^{l,r}(\omega) \quad. \end{equation}

For one ear (left) this can be interpreted as a frequency dependent ambisonic beamformer \begin{equation} s^l(\omega) = \chi_{nm}(\omega) [\breve H_{nm}^{l}(\omega)]^T \quad . \end{equation}